Critical Dynamics of a Two-dimensional Superfluid near a Non-Thermal Fixed Point 
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Critical dynamics of an ultracold Bose gas far from equilibrium is studied in two spatial dimen- 
sions. Superfluid turbulence is created by quenching the equilibrium state close to zero temperature. 
Instead of immediately re-thermalizing, the system approaches a meta-stable transient state, char- 
acterized as a non-thermal fixed point. A focus is set on the vortex density and vortex-antivortex 
correlations which characterize the evolution towards the non-thermal fixed point and the departure 
to final (quasi-)condensation. Two distinct power-law regimes in the vortex-density decay are found 
and discussed in terms of a vortex unbinding process and a kinetic description of vortex scattering. 
A possible relation to decaying turbulence in classical fluids is pointed out. By comparing the results 
to equilibrium studies of a two-dimensional Bose gas, an intuitive understanding of the location of 
the non-thermal fixed point in a reduced phase space is developed. 

PACS numbers: 03.75.Kk, 03.75.Lm 47.27.E-, 67.85.De 



I. INTRODUCTION 



Generically, the properties of an interacting many- 
body system far from equilibrium can change violently 
in time. This is not the case, however, when the system 
reaches a non-thermal fixed point. In the vicinity of such 
a fixed point the dynamics is expected to slow down and 
correlation functions to exhibit universal scaling behav- 
ior. Such transient states have been intensively studied in 
the context of classical turbulence [1-4] . Similar phenom- 
ena appear in degenerate quantum many-body systems, 
e.g., in superfluid helium or dilute ultracold quantum 
gases where superfluid or quantum turbulence (QT) has 
been discussed in great detail [5-12]. More general types 
of turbulence, so called wave turbulence have been stud- 
ied [13, 14], mainly in the framework of kinetic theory. 
Recent applications include far-from-equilibrium quan- 
tum systems, such as dilute ultracold Bose gases [15-25], 
the inflating and reheating early universe [26-30], and 
quark-gluon matter created in heavy-ion collisions [31- 
37]. Thereby, an extension of kinetic wave-turbulence by 
non-perturbative quantum-field-theory methods lead to 
the notion of a non-thermal fixed point [27, 28, 30, 33], 
in analogy to fixed points describing equilibrium as well 
as dynamical critical phenomena [38]. 

It was demonstrated in Refs. [16, 24] that in a 2- or 
3-dimensional superfluid Bose gas, such a non-thermal 
fixed point is realized by a state with a dilute random 
distribution of vortices or vortex lines, respectively, of 
both, positive and negative circulation. This gave the 
relation to QT [5-9, 14] as well as a link between weak 
wave turbulence described by kinetic theory and topo- 
logical excitations in non-linear wave systems. 
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FIG. 1: Dynamical evolution of a two-dimensional superfluid 
near a non-thermal fixed point (NTFP). The sketch shows 
the equilibration process after a quench, in the space of in- 
verse coherence length lc and inverse mean vortex-antivortex 
pair distance In- The 'Dynamical evolution' illustrates tra- 
jectories of decaying superfluid turbulence starting from the 
time where vortices appear, t = ty, approaching the NTFP 
and finally evolving towards equilibrium. The line marked as 
'Thermal states' qualitatively illustrates these quantities for 
thermal configurations [39-41], featuring a steady decrease 
of inverse coherence with inverse mean vortex-antivortex dis- 
tance and including a Berezinskii-Kosterlitz-Thouless (BKT) 
phase transition. An unbinding of vortices of opposite circu- 
lation characterises the approach to the NTFP before finally 
all vortices decay, lo — t 0, to establish equilibrium phase co- 
herence, here at a temperature below the BKT transition. 



In this paper, we study the non-equilibrium dynamics 
of a two-dimensional Bose gas evolving towards and away 
from a non-thermal fixed point (NTFP), the stationary 
properties of which were discussed in Refs. [15, 16, 24]. In 
this process the appearance and decay of vortex excita- 
tions play a crucial role. We monitor the vortex density 
during equilibration of the turbulent gas and reveal a 
bimodal scaling behavior in time. By following vortex- 
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antivortex correlations, we show that this phenomenon 
is directly related to a non-equilibrium vortex unbinding 
process. Ultimately, vortex excitations evolve into an al- 
most random distribution which constitutes the universal 
scaling at the NTFP. 
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In contrast to classical turbulence, the decay of super- 
fluid turbulence is typically accompanied by a build-up 
of coherence and quasi-condensation [17-22, 42, 43]. In 
Fig. 1, we sketch the projection of this process onto the 
space spanned by the coherence length Iq and the mean 
inter- vortex pair distance lr>- Note that Iq is defined here 
in terms of the participation ratio, cf. Eq. (7), which is 
a measure for the width of the normalized first-order co- 
herence function and is less sensitive to noise in the tails 
of s (1) 0). In this way, the dynamical evolution towards 
and away from a non-thermal fixed point can be com- 
pared to the properties of near-equilibrium states of a 
two-dimensional degenerate Bose gas [39-41, 44-46]. Ar- 
rows mark the direction of the flow and indicate that 
critical slowing down occurs near the NTFP. An unbind- 
ing of vortices of opposite circulation occurs during the 
approach of the fixed point before finally all vortices de- 
cay to establish full equilibrium phase coherence. 

Metastable multi-vortex states or, in one spatial di- 
mension, solitary waves, are also known to appear from 
strong fluctuations in the vicinity of the normal-fluid to 
supcrfluid transition. Crossing such a transition by vary- 
ing an equilibrium macroscopic parameter like temper- 
ature at a certain rate is well known to induce defect 
creation. Their number depends on the coherence length 
at the point where the parameter variation ceases to be 
adiabatic [47, 48]. Experiments with ultracold Bose gases 
following such Kibble-Zurek-type protocols [49, 50] as 
well as generating superfluid turbulence [10, 11] are pur- 
sued with increasing effort and could serve to discover 
and study systematically non-thermal fixed points. 

Turbulence has served, since the seminal work of Kol- 
mogorov [2-4] , as one of the first phenomena to develop 
renormalization-group (RG) techniques out of equilib- 
rium. The effectively local transport processes in mo- 
mentum, i.e., scale space, which are at the basis of turbu- 
lent cascades immediately suggest themselves for an RG 
analysis [51]. For two-dimensional ultracold gases, the 
dynamical evolution in the vicinity of the Berezinskii- 
Kosterlitz-Thouless (BKT) critical point [45, 46] was 
studied in Refs. [52], also in terms of a perturbativc 
RG analysis. For early work see [38] and references 
cited therein. A more general set of perturbative non- 
equilibrium RG equations for the one-dimensional sine- 
Gordon model near the Luttinger-liquid fixed point was 
derived in [53] , and the route to a non-perturbative analy- 
sis also for the strong-coupling regime is provided by out- 
of-equilibrium functional RG techniques [54-56] which 
will be followed in a forthcoming paper. 



FIG. 2: (Color online) Panels A — £ illustrate the initially oc- 
cupied momentum modes ki n i (shaded discs) with n(ki n i,£ = 
0) S> 1, for 5 different mean field configurations. The areas of 
the discs are proportional to the mean number of particles. A 
sixth initial condition A* is geometrically the same as A, but 
with initially occupied momenta k* ni = 4ki n i- If not stated 
otherwise, on average half of the total number of particles 
occupies the zero mode. 



II. DYNAMICAL SIMULATIONS 

In this paper we focus on the dynamical evolution of a 
two-dimensional dilute Bose gas towards a non-thermal 
fixed point and away from it to thermal equilibrium. We 
statistically simulate the far-from-equilibrium dynamics 
in the classical- wave limit of the underlying quantum field 
theory. The classical equation of motion for the complex 
scalar field </>(x, t) reads 



i<9 t ^(x, t) 



-^+. 9 |0(x,t)| 2 



(x,t). 



(1) 



Here, m is the boson mass, g quantifies the interaction 
strength in d = 2 dimensions, and, in our units, H = 1. 
Our computations are performed in a computational box 
of size L 2 on a grid with side length L = N s a s , lattice 
spacing a s , and periodic boundary conditions. We de- 
fine the dimcnsionless variables g — 2mg, t = t/r with 
lattice time unit r = 2ma 2 s and ip n (t) — V'na s exp(2ii), 
see [24] for further details. All simulations are per- 
formed with parameters g = 3 x 10 -5 and N/N 2 = 1525, 
where N is the total number of particles. If appropriate, 
we express length scales in units of the healing length 
£ = (2mgN/L 2 )^ 1 ^ 2 = 4.6 a s . We drop overbars in the 
following. We choose initial states with a few macroscop- 
ically occupied modes in momentum space, as illustrated 
in Fig. 2. Fluctuations around these mean values are 
introduced by sampling the initial field modes accord- 
ing to Gaussian Wigner distributions. Such statistical 
simulations, being also done under the name truncated 
Wigner approach, are quasi-exact in the classical wave 
regime of macroscopic occupation numbers. The type of 
initial conditions shown in Fig. 2 induces transport of 
particles and energy, which leads to vortex creation and 
turbulence. These field configurations describe quenched 
supcrfluid states, which can for instance be prepared by 
Bragg scattering of photons from a (quasi-)condensate. 
By varying the number and geometry of initially occu- 
pied modes, we can probe the initial-state dependence of 
our observablcs. 
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FIG. 3: (Color online) Left panel: Single-particle occupa- 
tion numbers n(k, t) = (|0(k,£)| 2 ) as a function of time t 
(in lattice units), for different discrete momentum modes 
k = (kn,0), k n = 2sin(n7r/JV s ) along the fc^-axis. Results 
are computed on a grid of size N s = 512, from initial field 
configuration A. Note the double-logarithmic scale. Lines 
show different power-law evolutions ~ t 2n . Right panel: Zero- 
mode occupation number n((ko, 0), t) and coherent popula- 
tion Ac = / d 2 a; (0(x, £)) | 2 as a function of time t (double- 
log, scale), for an average over 100 runs, grid size N s — 512 
and initial condition A. ty marks the time of vortex creation. 



A. Creation of vortices and turbulence 

In the presence of a non-vanishing coupling g the ini- 
tial states depicted in Fig. 2 are far from thermal equi- 
librium. During the first stages of the evolution coherent 
scattering into higher excited modes dominates. In Fig. 3 
(left panel) , we show the time evolution of the ensemble- 
averaged single-particle momentum occupation numbers 
n(k,t) = (\(j>(k, t)\ 2 ) as a function of time, for several 
momenta k = (ki,0) along the fc^-axis. One observes a 
power-law growth of momenta with k x > until t ~ 10 2 . 
This process, exhibiting fast power- law growth ~ t 2n of 
occupations can be understood from analytic mean-field 
calculations by approximating strong initial occupations 
to be time-independent. It is present for all our initial 
conditions and independent of spatial dimension. 

In Fig. 3 (right), we continue to follow the time evo- 
lution of the condensate mode, n(0,t). The decay of 
the zero-mode occupation is part of a nonlocal energy 
and particle transport to higher momenta. In coordinate 
space, this process leads to the formation of shock waves, 
which decay into large numbers of vortices, see Ref. [24] 
and videos of the evolution [57]. In addition, we study 
the coherent population Nq = / d 2 x |(</>(x, £))| for an 
initially phase coherent ensemble Nc(t = 0) = N. The 
dynamics preserve coherence until vortices form around 
time t = £y — 10 3 . This can be understood by consider- 
ing that the local phase angle y(x, t) of the complex field 
4> = |0|cxp{i</?} is determined by the positions of the vor- 
tices. Since vortices strongly interact, their trajectories 
in position space quickly randomise. Hence, in the en- 
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FIG. 4: (Color online) Vortex density p as a function of time 
t (in lattice units). Note the double logarithmic scale, (a) 
Evolution for various initial conditions as given in Fig. 2, 
averaged over 20 runs on a grid of size N a = 1024. Lines 
show different power-law evolutions. Assuming the area of 
a vortex to be given by Vv = (2£) 2 7r, a maximally dense 
packing of vortices would correspond to the vortex density 
Pmax = 1/VV — 1/(4tt£ 2 )- (b) Evolution for different initial 
zero-mode populations n(0)/N. Averages were taken over 20 
runs on a lattice of size N s = 1024, for initial condition A. 
Fitted parameters are c\ = 0.0026 and C2 = 0.0012. The clos- 
est approach to the NTFP is reached at t ~ (5 ... 10) X 10 5 . 

scmblc average the coherent population decays. Beyond 
this time, density fluctuations are significant only at mo- 
menta larger than the inverse healing length, k > l/£, 
whereas long-range fluctuations of the Bose gas are dom- 
inated by vortical flow. For t > ty, the zero-mode pop- 
ulation n(0, t) starts to increase, signaling the onset of 
phase ordering associated with vortex annihilations. This 
process will be studied in the following sections. 

B. Vortex density decay 

After the creation of vortices, the dynamical evolution 
exhibits a dual cascade in momentum space, transport- 
ing particles from intermediate to small momenta and 
energy from intermediate to large momenta [24]. The 
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single-particle momentum spectrum develops a quasi- 
stationary bimodal scaling, with characteristic exponents 
corresponding to the respective cascade processes. The 
system approaches a non-thermal fixed point [15, 24]. 
The low-momentum scaling of the single-particle momen- 
tum distribution can be related to the presence of ran- 
domly distributed vortices and antivortices [24]. In the 
present article, the evolution towards and away from the 
non-thermal fixed point is investigated. First, we will 
show that the approach to the non-thermal fixed point is 
accompanied by a change of the characteristic scaling of 
the ensemble averaged vortex density p(t) with time. 
Fig. 4a shows the time evolution of the vortex density 



(a) 



p(t) = (N v (t) + N A (t))/V 



(2) 



with iV v ( A ) (t) being the number of vortices (antivortices) 
in the volume V at time t, found in simulations start- 
ing from the initial conditions defined in Fig. 2. Vor- 
tices are counted by detecting their characteristic density 
and phase profiles. In all runs, vortex formation occurs 
around fry — 10 3 , apparent from the steep increase of 
vortex density around this time. For t > ty, two dis- 
tinct stages in the vortex density decay are observed, a 
rapid early stage and a slow late stage. Specifically, the 
vortex density follows power laws p(t) ~ t~ ai with two 
different exponents on, i = 1,2. The exponent during 
the early stage depends considerably on initial conditions 
1 < ol\ < 2, whereas the late stage features a decay expo- 
nent in a narrow interval 0.3 < a 2 < 0.4. From our anal- 
ysis given in Sect. II G below, we estimate that the closest 
approach to the NTFP is reached at t ~ (5 . . . 10) X 10 5 . 

We have repeated our simulations on various grid sizes, 
N s S {256, 512, 1024, 4096}. Thereby, we found that de- 
cay exponents saturate for and above N s = 512. We at- 
tribute deviations on smaller grids to effects from regular 
(intcgrable) dynamics of few- vortex systems [58] . We re- 
mark that the onset of the slow decay coincides with the 
development of a particular scaling behavior in the single- 
particle momentum distribution n{k) ~ fc -4 , which in 
Refs. [15, 16, 24] was shown to signal the approach to 
the non-thermal fixed point and the formation of a set 
of randomly distributed vortices. In this context, the re- 
duction of the vortex density decay exponent, compared 
to the early stage of rapid decay, is interpreted as a (crit- 
ical) slowing down of the nonlinear dynamics near the 
non-thermal fixed point. 

As shown in Fig. 4b, the vortex density decay at late 
times is not always given by a power law. By consider- 
ably increasing the initial population of the zero-mode, 
c. g. n(0)/N € {0.6,0.8,0.9}, we find that for some time 
the vortex density is better described by an inverse In 
function p(t) ~ l/ln(i). However, at a late times t > t* , 
the decay seems to converge to a power law from the 
slow-decay regime. By analyzing the dynamics of the 
single-particle momentum distribution, we could identify 
the time t* to be the time when compressible excita- 
tions have thermalized the high-momentum tail of the 
spectrum (for details see Fig. 10 in Ref. [24]). This is 
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FIG. 5: (Color online) Normalised vortex-antivortex correla- 
tion functions <?va defined in Eq. (3) as a function of radial 
coordinate r, for six different times U (in lattice units), (a) 
<?va(t") at times ti, i = 1,2,3 during the rapid-decay stage, 
averaged over 174 runs on a grid of size N a — 1024, using ini- 
tial condition A. Inset: Vortex density p as a function of time 
t, taken from the simulations for Fig. 4a. (b) gvA(r) at times 
ti, i — 4,5,6 during the slow-decay stage, averaged over 174 
runs on a grid of size N a = 1024, using initial condition A. 
Inset: Mean vortex-antivortex pair distance In as a function 
of time, from the simulations for Fig. 4a. 



consistent with the observation that the inverse-ln decay 
could not be observed for initial conditions with small 
zero-mode occupation where the high-momentum ther- 
malization happens faster. 



C. Vortex correlations 

In the following, the dynamical transition in the vortex 
annihilation dynamics is discussed in terms of character- 
istic features of the vortex-antivortex correlation function 



5VA(x,x',f) 



(p v (x,t)p A (x',t)) 



(3) 
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FIG. 6: (Color online) Sketch of different scattering events 
between two vortex- ant ivort ex pairs in d — 2 dimensions. 
Left panel: Scattering of two vortex pairs resulting in mu- 
tual annihilation of two vortices and the emission of density 
waves. Right panel: Scattering of two vortex pairs, leading 
to a change in the vortex anti-vortex distance Zd and pair ve- 
locity v. We remark, that once Zd ~ £ a vortex pair decays 
rapidly under the emission of density waves. 



where p v ' A '(x, t) = X)i^( x — x iM) is the distribution 
of vortices (antivortices) at time t in a single run. For 
sufficiently large ensembles, c/va is a function of r = |x — 
x'| only. 

In Fig. 5a, we show the evolution of <7va(?", t) as a func- 
tion of r for different times during the fast-decay stage. 
At early times, one finds a strong pairing peak near r = 0. 
This peak gets quickly reduced and a hole is 'burned' into 
the correlation function near the origin, see Fig. 5b. Fol- 
lowing the time evolution of the spatial vortex distribu- 
tion we observe that this involves qualitatively different 
processes: Mutual annihilations of closely positioned vor- 
tices and antivortices occur under the emission of sound 
waves. Further separated vortices can approach each 
other in different ways as illustrated in Fig. 6. The scat- 
tering of two pairs can directly lead to the annihilation 
of one pair under the emission of sound waves. We con- 
sider this to include events where the pair distance of 
one dipole reduces below a certain threshold, so that it 
looks like a density dip rather than a vortex-pair. This 
density dip can still interact with other vortices, but it 
will quickly vanish. Alternatively, the scattering reduces 
the vortex-antivortex separation within one pair while it 
increases it within the other, in accordance with the On- 
sager point- vortex model [44] . We refer to this character- 
istic change in gvA(r) as a vortex unbinding process. The 
scattering of a closely bound vortex off an isolated vortex 
is not shown, because it is included as a collision between 
a closely and a loosely bound pair. At around the time 
£3 < t < £4 the power-law exponent of the vortex density 
decay changes to about a third of its previous value, see 
the inset of Fig. 5a. Next, we compute the mean vortex- 
antivortex pair distance Zrj, by averaging over distances 
between each vortex and its nearest antivortex. In accor- 
dance to the previous discussion, Zd grows continuously, 
exhibiting two characteristic stages, see inset of Fig. 5b. 
At times t > 10 4 , lr>(t) approaches the power-law solu- 
tion Zd ~ P _1 ^ 2 , as expected for uncorrelated vortices. 
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FIG. 7: (Color online) Contributions to the total energy as 
functions of time t (in lattice units), averaged over 174 runs. 
We show the interaction energy E- m t, compressible energy E c , 
incompressible energy Ei, and quantum pressure E q , as de- 
fined in App. A, derived on a grid of size N s = 1024 and 
for initial condition A defined in Fig. 2. Note the double- 
logarithmic scale. See Sect. II D for a discussion of the power- 
law evolutions. 



D. Energy equilibration 

The main result of the previous section was the relation 
of different stages in the evolution of the vortex density 
p(t) to characteristic features in the vortex-antivortex 
correlation function gyA (r, t) . The vortex density decay, 
was shown to be accompanied by a dynamical vortex un- 
binding. This finding can be supplemented by consid- 
ering the evolution of different energies contained in the 
gas. In Refs. [5, 59] it was suggested to decompose the 
kinetic energy into an incompressible and a compress- 
ible component, which, for conciseness, we give details 
of in App. A. In this way, contributions from vortical 
excitations can be separated from other excitations such 
as sound waves. In Fig. 7, we show the time evolution 
of different energy components for initial condition A. 
One observes that the decay of the incompressible en- 
ergy during the early-time stage can be estimated to 
follow a power-law ~ i~ 125 and in the late-time stage 
~ t~ - 2 . This decay happens considerably slower than 
the vortex density decay in the early- (late-)time stage 
~ i _17 (^ £~ ' 3 ), discussed in Sect. II B. As a result, the 
incompressible energy per vortex grows as ~ i 45 at early 
times and as ~ t 0A at late times. Since the energy of a 
vortex pair increases with distance, this is in agreement 
with the phenomenon of increasing vortex-antivortex pair 
distance Zq- 

At late times, compressible and quantum-pressure en- 
ergy components develop into an equipartitioned state, 
also observed in decaying supcrfluid turbulence starting 
from a Taylor-Green vortex configuration [60] . 



E. Possible relation to classical turbulence 
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in d = 2 dimensions 

In the following, we discuss a similarity between our 
results and findings in classical fluid turbulence. Great 
efforts have been made to investigate freely decaying 
turbulence in two-dimensional classical fluids, see, e.g., 
Refs. [62-69]. In this context, special focus was set on 
the decay of the enstrophy f2(i), which is related to the 
vorticity u = V x v of the velocity field v(x, t) by 



A* 



A 



(4) 



It was found by different methods that the long-time de- 
cay of the enstrophy is given by a power-law Q(t) ~ i~ 7 , 
with 7 ~ 0.35 - 0.4 [62-69]. 

In supcrfluids, vorticity is concentrated in the vortex 
cores. The vorticity of a turbulent flow consisting of M 
vortices with circulations Ki and positions Xj. for i < M, 
is given by tj(x.,t) = Yli K i^( x ~ X «W)- Hence, in a 
flow consisting of vortices with circulations Ki = ±1 the 
enstrophy reads 



f 



fi(t) = -5(0)M(i) , 



(5) 



which is, with 5(0) ~ 1/V, proportional to the vortex 
density, Q(t) cc p(t). As shown in Fig. 4a, in the late- 
time stage, our results are in accordance with the results 
from classical turbulence. However, we point out that the 
mechanisms of enstrophy decay in the two systems are 
fundamentally different. Whereas in superfluids vortices 
annihilate, the main process of vorticity decay in classical 
two-dimensional fluids is the merging into larger vortices. 

We finally remark that in numerical simulations of 
freely decaying classical turbulence a crossover between 
two stages of power-law decay similar to our findings has 
been reported in Ref. [69]. 



F. Kinetic theory of vortex scattering 

The decay of the vortex density has been investigated 
in two-dimensional classical fluids [62-69] and supcrflu- 
ids [43, 70-72], mainly in the presence of driving or dis- 
sipation. Several authors have proposed kinetic theo- 
ries building upon assumptions about the decay process 
[61, 62, 67, 68]. The vortex decay cannot be explained by 
a simple model of independent vortices and antivortices 
moving towards each other to minimize the energy. Ne- 
glecting interactions with sound vortex-antivortex pairs 
perform a collective motion perpendicular to their rel- 
ative distance vector without changing their distance. 
This motion quickly leads to pair-pair collisions and es- 
tablishes a kinetic-theory picture for vortex pairs. In 
Fig. 6, we showed two examples of vortex-pair scatter- 
ing processes altering the vortex density p(t) and vortex 
correlation functions gvA( r >0- 

Assuming that the vortices are moving in pairs and 
that annihilations happen as the result of collisions of 
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FIG. 8: (Color online) Vortex and antivortex positions shortly 
after turbulence creation at t ~ ty, for two single runs of the 
simulations on a grid of size N s = 1024. Left panel: Initial 
condition .4*. Right panel: Initial condition A, as defined in 
Fig. 2. 



vortex pairs, the decay rate for the number Ny of vor- 
tices follows from the number of dipoles iVn ~ Ny and 
dtN-£>(t) ~ — N-£,(t)/r, with mean free collision time r. 
The mean free collison time r = l m { p /v is given by the 
mean velocity of the pairs v and the mean free path 
^mfp = V/ (ct-/Vd) with cross section a. Both, mean veloc- 
ity and cross section depend on the number of vortices 
via the mean vortex-antivortex pair distance Id accord- 
ing to v ~ 1/Id and a ~ Zd- These considerations result 
in a rate equation 



d t Ny(t) = -cN v /t 



(6) 



with dimcnsionless constant c. Eq. (6) has the power-law 
solution Ny(t) ~ t~ l . This decay law can only be ob- 
served for specific initial conditions A* in the early-time 
stage. We attribute deviations from i -1 scaling during 
this period to an inhomogeneous distribution of vortices, 
encountered for certain initial conditions. To give an ex- 
ample, we show the vortex distributions created from the 
two different initial conditions A* and A in Fig. 8. 

During the late-time stage, only a few vortices are 
bound in pairs, while most vortices are loosely bound and 
interact equally with a larger number of vortices around 
them. We heuristically take this into account by consid- 
ering a modified scattering cross section a ~ luN^,. The 
resulting kinetic equation reads d t Ny(t) ~ —N y , with 
solution Ny(t) ~ i -1 / 3 , and hence yields the reduction 
of the decay exponent observed at late times in our sim- 
ulations. 



G. Phase correlations 

In the remainder of this article, we focus on the growth 
of long-range coherence at late times, associated with the 
annihilation of topological defects [17-22, 42, 43]. From 
this point of view freely decaying superfluid turbulence 
is a particular example of phase-ordering dynamics after 
a quench into the ordered phase [73]. Whereas in three 
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FIG. 9: (Color online) Trajectories of multi-vortex states in 
the space of inverse coherence length 1/lc and inverse mean 
vortex- ant ivort ex distance 1/Id, starting from t = ty. Arrows 
are added to guide the eye along the time direction. Dashed 
lines mark the minimal values 2/L = 0.018£ _1 , available on 
a grid of size N, = 512. Our understanding of the NTFP 
as a configuration with a few, maximally separated pairs on 
an otherwise maximally coherent background implies it to 
be located near the crossing of the dashed lines. Hence, the 
NTFP is approached most closely between t ~ 5 x 10 and 
t ~ 10 6 = Pi. Closed (red) circles mark an average over 174 
runs, N a = 1024, initial time t\ = 2.3 x 10 3 and final time 
tf = 1.3 x 10 , closed squares an average over 1223 runs, 
with N s = 512, ty = 1-7 x 10 3 and t t = 2.6 x 10 s ; open 
squares an average over 16 runs, N s — 512, ij = 2.6 x 10 , and 
tf = 4.2 x 10 6 . Note that the symbols are equally spaced on a 
logarithmic time scale. We indicate the times to — 1.7 x 10 3 , 
ti = 1.6 x 10 4 , ta = 1.3 x 10 5 , i 3 = 1.0 x 10 6 , U = 4.2 x 10 6 . 
The average number Ny(t) of vortices left in the system is 
N Y (t ) = 97.7, N v (ti) = 21.8, N v (t 2 ) = 11.7, Ny(t 3 ) = 5.8 
and N v (t 4 ) = 2.1. 



dimensions a second-order phase transition connects a 
normal- fluid and a superfluid phase, a Bose gas in two 
dimensions experiences a Berezinskii-Kosterlitz-Thouless 
(BKT) transition [45, 46]. For the two-dimensional ultra- 
cold Bose gas, experimental and theoretical results sup- 
port the understanding of the phase transition in terms 
of vortices undergoing an unbinding-binding transition 
[40, 41, 50, 74-77]. 

In this context, we are interested in a comparison 
between correlation properties observed in the non- 
equilibrium dynamics near a non-thermal fixed point 
and those known from equilibrium studies. We compute 
the dynamical trajectory of the vortex gas in the space 
of inverse coherence length and inverse mean vortex- 
antivortex pair distance. We compare our results to simu- 
lations of a thermal two-dimensional Bose gas specifically 
for our system parameters. For this, we evolve field con- 
figurations in time which are initially close to a thermal 
Rayleigh- Jeans distribution at a temperature T. After 
equilibration is reached, we compute the position of the 
states in the above phase space for different T. 

We define a coherence length Iq in terms of the 
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FIG. 10: (Color online) Multi-vortex states in the space of 
inverse coherence length and inverse mean vortex-antivortex 
distance, (a) Thermal configurations (Zr 3 1 (T),?Q (T)) for a 
range of temperatures T, increasing from bottom right to 
top left. The solid line marks the point where the decay 
of the g' 1 '(r)-function changes from algebraic to exponen- 
tial, signaling the BKT-transition. (b) Comparison of the 
thermal line (ip 1 (T), ij _1 (r)) for the same range of temper- 
atures T with the corresponding dynamical evolution. (Same 
data as in Fig. 9.) Dashed lines mark the minimal values 
2/L = 0.018f _1 , available on a grid of size N s = 512. Note 
that the (l/Zo)-axis interval [0.25,0.55] has been cut out. 

participation ratio [78] of the angle-averaged first- 
order coherence function g^(r) = j d9 (0*(x)<^>(x + 

r))/VMx))Mx + r)), 

lc=Ujdr\gW(r)A , (7) 

with Af = [JdrgW(r)]- 2 . It measures the spatial ex- 
tension of the first order coherence function. Other than 
r coh = / dr r 2 (r) / J drrg^'(r), the quantity lc does 
not sum up values of g^(r) weighted by the distance, 
which would enlarge insignificant contributions at large 
r. In addition, it gives meaningful results also in the 
case of large coherence g^'(r) ~ 1, where for instance 
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FIG. 11: (Color online) First order coherence function g^ 1 ' 
as a function of radial coordinate r. We show one example 
for the system during the equilibration process at time t = 
1.05 x 10 6 , with l c = 41.6$, l* c = 26.2$, and one for a thermal 
configuration, giving lc = 51.3$, Iq — 14.8$. 



the FWHM-measure can not be applied any more. Note 
that in equilibrium this quantity smoothly interpolates 
between the regime of exponential decay of above the 
BKT-transition, where the exponential coherence length 
$c is defined as g^'(r) ~ exp(— r/$c), and its power-law 
decay in the supcrfluid regime. 

In Fig. 9, we follow the time evolution of the gas for 
t > ty. One can observe, that a state of low coherence 
and small mean vortex-antivortex pair distance evolves 
towards larger coherence and larger vortex-antivortex 
separation. As discussed in Sects. II B and II C, this 
is due to vortex annihilations and vortex-antivortex un- 
binding. For times t > 10 4 , the coherence length grows as 
l c ~ p" 1 / 2 , in the same way as Zd shown in Fig. 5b. The 
evolution considerably slows down for 1/lc ~ — ^ 0. 
In this regime, the Bose gas shows characteristic scaling 
properties, see Refs. [16, 24], which indicate the pres- 
ence of the non-thermal fixed point [15]. After spending 
a long time near this point, the mean vortex-antivortex 
pair distance declines. This is a sign that the last remain- 
ing vortex-antivortex pairs reduce their size prior to their 
annihilation and the equilibration of the system. Around 
the same time the power-law in the vortex density decay 
shown in Fig. 4 breaks down. 

Our understanding of the NTFP as a configuration 
with a few, maximally separated pairs on an otherwise 
maximally coherent background implies it to be located 
near the crossing of the dashed lines. Hence, the NTFP 
is approached most closely between t ~ 5 x 10 5 and 
t ~ 10 6 = t 3 . 

To set the above evolution in relation to equilibrium 
configurations, we show, in Fig. 10a, the thermal line 
(Zp (T), Zq 1 (T)) for a range of temperatures T for which 
the system shows a non- vanishing zero-mode population. 
Note that, in order to define what counts as a bound 
vortex pair, we filter out field fluctuations on scales 
smaller than 0.55$ before detecting vortices and antivor- 
ticcs. Hence, the resulting inverse of the mean vortex- 



antivortex pair distance represents a lower bound, and 
the separation of the NTFP from the thermal configura- 
tions becomes obvious. 

In view of the thermal results it is useful to consider 
an alternative definition for the coherence length. The 
length Iq = J drg^'fr) shares the above mentioned ad- 
vantages of the participation ratio. In addition, it does 
not overestimate the coherence of flat distributions. The 
thermal g^'(r) functions show a fast decay at short dis- 
tances r s to a value and are almost constant for 
r > r s . Fig. 11 shows two typical gW-functions, one 
for the system during the equilibration process at time 
t = 1.05 x 10 6 , with l c = 41.6$, l* c = 26.2$, and one for 
a thermal configuration, giving lc = 51.3$, l c = 14.8$. 
Due to the small area under g^(r) up to r s the nor- 
malization M in the definition of the participation ratio 
enlarges the integrand of Eq. (7) close to unity. Hence, 
lc becomes large while l c is less affected by this effect. In 
Fig. 10b we compare the thermal line (l^ 1 (T) , l^" 1 (T)) 
for the same range of temperatures T with the corre- 
sponding dynamical evolution. 

In Fig. 1, our findings arc summarized qualitatively in 
a reduced phase space of the vortex gas. In this way, 
the 'Dynamical evolution' of decaying superfluid turbu- 
lence is compared to an estimate of the expected equilib- 
rium configurations illustrated by the dashed line marked 
'Thermal states' along which the temperature T varies. 
The most significant difference between these two lines 
is that the non-equilibrium dynamics is characterised by 
an increase of the mean vortex-antivortex pair distance 
with increasing coherence, whereas equilibrium configu- 
rations are expected to feature a decrease in pair dis- 
tance with increasing coherence. A slow-down of the dy- 
namics together with characteristic scaling of the single- 
particle momentum distribution, observed in the regime 
of large coherence and large vortex-antivortex pair dis- 
tance, marks the position of the non-thermal fixed point . 
Finally, when all vortices have annihilated, the system 
reaches the 'Thermal states'-line deep in the superfluid 
regime. 

III. CONCLUSIONS 

We have studied the non-equilibrium dynamics of a 
two-dimensional dilute ultracold Bose gas towards and 
away from a non-thermal fixed point (NTFP). Following 
an initial quench, evolution towards a fixed point ap- 
pears to be a generic feature of the (quasi-)condensation 
process and the build-up of coherence. In the course of 
a critically slowed evolution, vortex excitations evolved 
into an almost random distribution reflected in the scal- 
ing of the single-particle spectra at the NTFP. We showed 
that the vortex-density decay is directly related to a non- 
equilibrium vortex unbinding process. 

Our results allow to draw a picture of the NTFP. The 
evolution path towards and away from the fixed point is 
summarized schematically in Fig. 1. There, we sketch the 
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evolution of the multi- vortex states in the plane spanned 
by the coherence length Iq, defined in terms of the partic- 
ipation ratio, and the mean inter- vortex pair distance Id. 
This picture allows to compare our results with quasi- 
equilibrium studies of a two-dimensional vortex gas. The 
NTFP emerges to bear similarities to the equilibrium 
Berezinskii-Kosterlitz-Thouless (BKT) fixed point. The 
NTFP is characterized by a few, in the extreme case, one 
pair of far separated anti-circulating vortices. While the 
BKT transition also features unbinding of vortices, finite 
temperature implies the simultaneous excitation of many 
rotons, i.e., strongly bound vortex-antivortex pairs. The 
NTFP is identified by strong wave turbulent scaling in 
the infrared limit [15], n(k) ~ k~ 4 . At the same time, the 
high-energy modes can be populated in a much weaker 
way, e.g., at considerably lower temperature than the 
BKT critical temperature or remain out of equilibrium. 
The details of the UV mode populations are determined 
by the way how the NTFP is being approached. We 
emphasize that the approach of the NTFP is a generic 
but out-of-equilibrium process and that eventually, the 
system will decay to an equilibrium state potentially far 
away from the NTFP. In the examples discussed in this 
paper, the total energy is sufficiently low such that the 
final equilibrium state emerges to be considerably below 
the BKT critical temperature. 

The way we force the system here as in the work re- 
ported in [16, 24] to approach the NTFP generalizes that 
of Kibble and Zurck. A strong sudden quench replaces 
the more or less adiabatic approach of the BKT phase 
transition. To what extent the BKT phase transition 
can be understood as happening within a class of ther- 
mal states near the NTFP studied here needs to be clari- 
fied by analyzing full out-of-cquilibrium rcnormalisation- 
group equations in comparison with standard descrip- 
tions in thermal equilibrium. We point out that the con- 
cept of an NTFP is far more general than the specific sit- 
uation studied here, see, e.g., Refs. [15, 25, 27-30, 32, 33]. 

Further important questions concern the relation to 
fully developed quantum turbulence [5-12] which is be- 
lieved to exhibit a quasi-classical Kolmogorov-Obukhov 
scaling [2, 3] of the radial energy, E(k) ~ fc -5 / 3 , in the 
infrared regime below the mean inverse distance between 
vortices [5-9]. 

In a recent experiment, a strong reduction of the vortex 
decay rate at late times has been reported [12]. Quanti- 
tative experimental observation of the predictions made 
here could give new insight into the character of the 
NTFP and its relation to quantum turbulence. 
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Appendix A: Quantum hydrodynamic energy 
decomposition 

In this appendix we give details of the energy decom- 
position introduced in Refs. [5, 59]. To exhibit vorti- 
cal flow and define the decomposition we use the po- 
lar representation </>(x, i) = y / n(x, t) exp{i<^(x, t)} of the 
field in terms of the density n(x, £) and a phase an- 
gle <p(x.,t). This allows to express the particle current 
j = i((f)*\7(f> — (f>\7<j>*)/2 = nv in terms of the velocity 
field v = V(p. With this, we decompose the kinetic- 
energy spectrum following Ref. [5, 59], splitting the to- 
tal kinetic energy £ kin = J d d x (|V(/>(x, i)| 2 )/(2m) as 
i?kin = E v + E q into a 'classical' part and a 'quantum- 
pressure' component, 

E v = ^- fd d x(\^\ 2 ), (Al) 
Ira J 

Eq =^LJ d Mi v ^i 2 >- ( A2 ) 

The radial energy spectra for these fractions involve the 
Fourier transform of the generalised velocities w v = \/ri\ 
and w q = V-y/n, 

E s (k) = ^fk d - 1 dn d (\ Ws (k)\ 2 ), S = v,q. (A3) 

Following Ref. [5, 59], the velocity w v , is furthermore 
decomposed into 'incompressible' (divergence free) and 
'compressible' (solenoidal) parts, w v = w; + w c , with 
V ■ w; = 0, V x w c = 0, to distinguish vortical supcrfluid 
and rotationlcss motion of the fluid. 
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